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Abstract. We analyze numerically a forward-backward diffusion equation with a cubic- 
like diffusion function, -emerging in the framework of phase transitions modeling- and its 
"entropy" formulation determined by considering it as the singular limit of a third-order 
pseudo-parabolic equation. Precisely, we propose schemes for both the second and the 
third order equations, we discuss the analytical properties of their semi-discrete counter- 
parts and we compare the numerical results in the case of initial data of Riemann type, 
showing strengths and flaws of the two approaches, the main emphasis being given to the 
propagation of transition interfaces. 

Keywords. Phase transition, pseudo-parabolic equations, numerical approximation, 
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1. Introduction 

The aim of the present article is to investigate numerically the solutions to a nonlinear diffusion 
equation of the form 

du d 2 d>(u) 

(1) m=d^ (x,0eRx(0,+oo) 

where u : R x [0, +oo) — > R, for a non-monotone diffusion function <p : R — > R; specifically, the 
function cf> is to be cubic-like shaped, i.e. monotone decreasing in a given interval / and monotone 
increasing elsewhere. 

Due to the presence of intervals where (p is a decreasing function, the initial value problem 
for ([T} is ill-posed and an appropriate (non classical) definition of solution has to be introduced. 
In this respect, a common paradigm states that the loss of a well-posed classical framework may 
emerge as consequence of some simplification of a given complete physical model, consisting in 
neglecting some higher order term because of the presence of some small multiplicative factor 
£. The usual strategy is then to consider classical solutions to the original higher order equation, 
pass to the limit as a definite parameter s tends to zero and regard such singular limit as the 
"physical" solution of the problem under consideration. Whenever possible, it is interesting to 
prescribe a definition of generalized solution to the reduced problem that is consistent with the 
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limiting procedure and, at the same time, is determined by proper conditions to be verified by such 
a generalized solution, without explicit reference to the disregarded higher order terms. In this 
respect, a classical example is given by scalar conservation laws where solutions can be defined 
through a vanishing-viscosity approach and well-posedness is restored by considering an entropy 
setting. Different higher perturbations may, in principle, lead to different limiting formulation, 
as in the case of van der Waals fluids, where viscosity and viscosity-capillarity criteria give raise 
to different jump conditions in the limit. E71 . In the same spirit, the approach through singular 
perturbations is commonly used in the context of calculus of variations to determine appropriate 
minima of non-convex energy functionals (among others, see the classical reference Ref. [16]). 
A related approach consists in studying the structure of the global attractors of the higher order 
equations in the singular limit, as it has been performed for a relaxation viscous Cahn-Hilliard 
equation in Ref. |9). 

Here, the model ([T]) arises as a reduced equation in the context of phase transitions and we 
consider the formulation for ([T]), established in Refs. EH . 112211 . 151 . obtained as singular limit of 
the third order equation of pseudo-parabolic type 



originally considered in Ref. ifTTl . The denomination "pseudo-parabolic", relative to the presence 
of the time derivative of the elliptic operator d /dx , has been proposed in Refs. Eoll . Il29lk equa- 
tions of the same form are also called of Sobolev type or Sobolev-Galpern type. A sketch of the 



order fluids (see, respectively, Ref. [3], Ref. [4] and descendants). Besides, the third order term in 
([2]) appears in the description of long waves in shallow water (see Refs. |2j, [20] and descendants), 
as an alternative to the usual purely spatial third order term considered in the Korteweg-DeVries 
equation. Taking the singular limit of such higher order equations gives in general a different 
formulation with respect to the classical entropy formulation obtained by means of the vanishing 
viscosity limit,[31]. Regardless of the physical context, the use of the third order term as a reg- 
ularizing term for an ill-posed diffusion equation of the form ([T]), is called the quasi-reversibility 
method, and it has been used by many authors, mainly in the linear case (see Ref. Ifl9l and refer- 
ences therein). 

Independently of the monotonicity of the function cf>, the Cauchy problem for equation Q is 
well-posed, IfTTl : moreover, the singular limit u of the family {u E }, solutions to ([2]), can be described 
in term of Young measures, [21 ] [22]. The sign of <p' is relevant for the dynamical properties of 
the solutions: precisely, the constant state u turns to be (asymptotically) stable if (p'{u) > and 
unstable if (p'{u) < 0. In the case under consideration, it is possible to distinguish two (unbounded) 
stable phases divided by a (bounded) unstable phase, the so-called spinodal region. Generically, 



(2) 





derivation of such an equation in the phase transition setting is given in Section [2] 

Equations similar to (|2]) arise also in modeling fluid flow in heat conduction and shear in second 
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initial data with values in the spinodal region generate oscillations which lead to patterns mixing 
the two stable phases. 

Assuming additionally that the limit u of the family {u E } is piecewise smooth and that it takes 
values only in the regions where the function (p is monotone increasing, it is possible to determine 
appropriate admissibility conditions, called entropy conditions, on a function in order for it to 
be limit of solutions to the higher order equation Q (see Proposition 2.1 below). For such class 
of generalized solutions, from now on denominated two-phase solutions, uniqueness and (local) 
existence have been proved, |[T4l . 

Here our aim is to investigate numerically the behavior of solutions to ([T]) exploring two pos- 
sible strategies: on the one hand, by discretizing the third-order equation Q and considering the 
behavior of the algorithm as e — > 0; on the other hand, by providing an approximation of the solu- 
tions to Q, taking advantage of the entropy conditions. Schematically, if we denote by h the size 
of the discrete spatial mesh and by s the parameter appearing in Q, both values are assumed to 
be small, ideally tending to zero: in the first approach, we let e — > + and then h — > + . Viceversa, 
in the latter, we take the limits in the opposite order. Absence of uniformity of one parameter with 
respect to the other implies that the two procedures lead, in principle, to different results and we 
are interested in comparing them. Specifically, dealing directly with (|2]) obliges to treat solutions 
taking values also in the region where (j> is decreasing, with the consequent appearance of strong 
oscillations due to the exponential instability of constant states for the backward heat equation. On 
the contrary, the two-phase solutions are designed to take values only in the region of stability {<p 
increasing) and to cross the unstable region in single specific transition points where appropriate 
entropy transmission conditions should be satisfied. 

We benchmark the two approaches against Riemann problems, determined by the values (u~, u + ) 
with belonging to the two different stable phases. Depending on the location of such values, 
the initial phase transition discontinuity, hereafter called interface in reference to the porous media 
equation [15], may or may not move (see Section|2]). As expected, relevant differences appear in 
the location of phase boundary depending on the algorithm chosen (see Figure[T] and Section[5]). 

Coming back to the original modeling, it would be more appropriate to consider equation ([T]) 
as the singular limit of some higher order equation with a much more complicated structure with 
respect to the one of ((2]) (see Section [2]). To our knowledge, in this case there is still no available 
limiting formulation, meaning that one does not know which are the admissibility conditions for 
the solutions to ([T]) to be singular limits of solutions to the corresponding higher order equation. 
Hence, the only possible approach consists in taking the singular limit after having discretized the 
complete equation. Such procedure can not be rigorously justified because of the non-uniform 
dependence of s with respect to h. Therefore, analyzing a simplified situation where the two 
approaches are both available and can be compared, can be heuristically considered as a reliability 
test on the exchange of the limits e — > + and h — > + . 



Figure 1 . Evolution of the relative error of the position of the interface for an initial data 
of Riemann type with (u~, u + ) = (-2,4) for the two approaches: explicit scheme for the 
pseudo-parabolic equation (blue), two phase scheme (red). 

The article is organized as follows. In Section |2| we sketch the derivation of the modeling 
equation in order to identify the regime of parameters that leads to ([2]) and ([T]). Also, we recall 
the main results on the analytical properties of the equation, with particular care to the two-phase 
formulation. We also deduce an explicit formula for the solution of the Riemann problem in the 
case of a piecewise linear diffusion function cf>. In Section[3j we propose a standard finite-difference 
algorithm for the third order pseudo-parabolic equation ((2]), as considered in Ref. J5), and we 
analyze the qualitative behavior of the solutions for the corresponding semi-discrete scheme. In 
Section[4j we test the fully-discretized algorithm in the case of Riemann type initial data, showing 
the emergence of incorrect oscillations in the solution. Finally, in Section [5] we implement an 
algorithm for the two-phase solution of ([T]), consisting in coupling finite-differences schemes for 
the two stable phases with some appropriate transmission conditions at the phase boundary and 
we compare the results given by the two approaches. 

2. Physical and analytical background 

Phase transitions modeling is a wide area of active research in pure and applied mathematics. 
The basic target is to describe the behavior of a mixture with many different (stable or unsta- 
ble) phases, with special attention to the dynamics of separating interfaces and pattern formation. 
Various approaches and models have been proposed and analyzed, depending on the specific con- 
text. Probably, the most celebrated is the Cahn-Hilliard model, which is based on defining an 
appropriate cost functional which determines two preferred states and, at the same time, penal- 
izes inhomogeneities by means of a gradient term (for a survey on Cahn-Hilliard and phase field 
models, see Ref. lloll). 
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Here, we consider a different model, devised in Ref. ifTTl . that still preserves the basic structure 
of two stable "competing" phases and gain regularity from memory effects, without taking into ac- 
count the cost of transition layers. In this Section, we present a simplified version of the derivation 
of the model equation, and then we recall the main known analytical results. 



Derivation of the model. Consider a generalized diffusion equation, in a one-dimensional domain 
described by the space variable ieR, 

du d 2 v m 

(3) - = 5? «R.,>0. 

where u = u(x,t) and v - v(x, t) are real valued functions. Following Ref. 021 > the unknown v, 
called potential, is supposed to satisfy the relation 

(4) v := tf/(u) + f 9'(t - s) (if/(u) - <p(u))(x, s) ds 

J — CO 

or, equivalently, 

C d 

v = (p(u) + 9{t - s) — (4>(u) - <p{u)){x, s) ds, 

J-oo OS 

where 9 is a relaxation function, assumed to be monotone decreasing and such that 9(0) = 1 
and 9(+oo) = 0. The integral term in the definition of v represents memory effects localized in 
space, with decay described by the relaxation function 9. The case with no-memory is obtained by 
choosing -9' equal to the usual Dirac distribution concentrated at zero. In this case, there holds 
v = <p(u) and equation ([3]) reduces to a standard nonlinear diffusion equation. 

The relaxing structure Q mimes the analogous expression used as a stress-strain law in the 
modeling of viscoelastic media and based on the Boltzmann supeiposition principle (for the gen- 
eral mathematical theory, see Ref. [25] for viscoelastic materials with memory and Ref. [24] for 
viscoelastic flows). In that context, functions if/ and (p represent, respectively, the instantaneous 
elastic and the equilibrium stress-strain responses. Similar interpretations can be given also in the 
present context as a contribution to the potential v. 

The simplest choice for 9 = 9(t) is a single exponential function e~^ T , t > 0. In this case, it 
is possible to write a partial differential equation for the unknown u. Indeed, since 9(0) = 1 and 
9' - -9 It, there holds 

dv dilf(u) 1 C 1 d . , M(u) 1 

9(t - s) —(i//(u) - (f>(u))(x, s) ds - — (v - (p(u)) 



-f 

T J-o 



dt dt t J.co ds dt t 

Therefore, a relaxation function of exponential type corresponds to the assumption that the poten- 
tial v solves a relaxation-type equation 

d 1 

- (v - i/f(u)) - - (0(11) - v) 
dt t 
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The case with no-memory effects is obtained by taking the singular limit t — > and, again, 
formally setting v = <p(u). Substituting in Q differentiated with respect to t, we infer 

d 2 u d 2 (dip(u) 1 \ 1 d 2 v 



dt 2 ~ dx 2 \ dt + r m J rdx 2 ' 
hence, using ([3]) to eliminate the second derivative of v, we obtain 

d 2 u du d 2 I dtl/(u) 

T—r + — = —r\(f>(u) + T 



dt 2 dt dx 2 \ rw dt 
Incorporating in the analysis a van der Waals term to take into account the free-energy cost of in- 
homogeneities would lead to a relaxation viscous perturbation of the classical fourth-order Cahn- 
Hilliard equation (containing also a fifth-order term arising from the coupling of the inhomogene- 
ity cost term and the memory effect). 

From now on, we choose tp in the form ip(u) - ku with k > 0. Hence, the unknown u satisfies 
the equation 

d 2 u du d 2 I du 

(5) T w + di=M m+£ * 

where e := tk. In Ref. [17], the regime t — » + with e being fixed has been considered, yielding 
the pseudo-parabolic equation Q. The nonlinear diffusion equation ([T]) arises in the limiting 
regime s — > + of equation ((2]). Equivalently, we can consider the system of partial differential 
equations 

( du dw 
dt dx 
dw d(p{u) d 2 w 
dt dx dx 2 

where w denotes the space-derivative — of the potential v, and then consider equations pb and 

dx 

([T]) as formally obtained by passing to the limit (in this order) t — > + and e —> + . 

Let us stress that the assumptions on the relative sizes of parameters r and s considered is 
not motivated by any physical argument and, thus, is a pure mathematical simplification. To our 
knowledge, a complete description of the limiting behavior of the equation (|5]> as t, e — > + for 
different relative sizes of the two parameters is still not available. 

Overview on the analytical theory. Next, let us review the main analytical results known about 
([T]) and ([2]) and their relations. Let cp : R — > R be a locally Lipschitz continuous function with a 
cubic-like graph with a local minimum with value A and a local maximum with value B, for some 
A < B: precisely, we assume 

(p strictly increasing in (-oo, b) U {a, +oo), (p strictly decreasing in (b, a) 

for some b < a with <p(a) - A and cp{b) - B. For later reference, we assume the existence of 
c e (-oo, b) and d e (a, +oo) such that 0(c) - A and <p{d) - B. 



7 



The monotonicity of the function cf> distinguishes three different phases 

(6) phase S~ = (-00, b], phased = (b, a), phase S + = [a, 00). 

Being interested in the dynamics generated by the different stability properties of the regions 
(-00, b), (b, a) and (a, +00), we will restrict the analysis to the case of a piecewise linear function 
and, specifically, to the case 

3 

(7) 0(h) = 2b+-(|1-h|-|1+b|). 

In this case c = -2,b = -I, a = l,d = 2 and A — —l,B - 1 (see Figure [2]). The third-order 




Figure 2. The piecewise linear diffusion function if>. 

term in ((2]) acts as a regularizing term and the initial value problem on a bounded interval /, with 
Neumann boundary conditions, for the pseudo-parabolic equation ([2]) is well posed in the classical 
sense for initial data either bounded or continuous, iTTTIl lfl"8l . In Ref. [ 17], it is also shown that 
any bounded measurable function u = u(x) such that <p(u) is constant is a stationary solution of the 
equation and, if, in addition, <p'(u(x)) > for any x, such a solution is also asymptotically stable 
(with exponential decay rate) with respect to zero-mass initial perturbations. In particular, given 
u~ < b < a < u + and xq € /, step functions 

(u~ x < Xq, 

+ 
U X > Xq, 

satisfying the transmission condition 

cp{u~) = (p{u + ) e (A,B) 

are attracting time-independent solutions to ([2]). In the limit s —> + such functions persist to be 
solution of the corresponding limiting equation ([T]). 

In Refs. ETI . E2ll . a deep analysis of the limit e — > + of the solutions u E to the pseudo- 
parabolic equation (|2]) has been performed within the Young measure framework. Poorly speak- 
ing, the analysis of Plotnikov shows that the singular limit of the family of solutions u E can be 
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represented as 

(9) u = A-(<f>-)-\v) + A°(<p°)- l (v) + A + (cf> + )-\v) 

where A , A* and v are bounded functions and <f>~, cfp, + represent the restriction of <p in (-00, b), 
[b, a] and [a, +00), respectively. Functions A , A ± are non-negative and satisfy 

[ A~(x,t) - 1 if v(x,t) < A, 



A (x, t) + A ( \x, t) + A + (x, = 1, and 



A + (x,t) = l if v(x,t)>B, 



for any (x, t) under consideration. Finally, for any non decreasing function g e C l (R), setting 
G* := A~ G((f)" 1 (v)) + A°G((4> )-\v)) + A + G(((p + )- l (v)), where 

G(s):= f g(<P(o-))dcr, 
Jo 

the following inequality holds in the sense of distribution 

dG* d I dv\ , (dv\ 2 

Representation (|9]> can be heuristically interpreted by stating that, in the limit e — > 0, the solution 
u is decomposed in the sum of three contributions, each relative to one of the three phases S ± and 
H. Functions A , /I* represent the fraction of u in the corresponding phase. Inequality ( fTO] ) can 



be interpreted as an entropy inequality, recording the irreversibility of the limiting process. At the 
present time, it is not known if all these conditions are sufficient for determining uniqueness of the 
Cauchy problem (results on large-time behavior are obtained in Ref. |28 ]). Hence, the previous 
list of conditions is unsatisfactory for determining a well-posed limiting dynamics. 

Two-phase solutions. Linearization around an unstable value u € (b, a) readily shows that solu- 
tions of the third order equation Q tend to escape from the unstable interval exponentially fast, 
precisely with rate exp(-cp'(u) tie). In the limit s — > + , values in (b, a) instantaneously exit from 
the unstable region; hence, solutions to ([T]) are expected to take values generically in the stable 
regions. Thus, a natural point of view is to consider solutions composed of "pure" stable phases, 
meaning that, in the decomposition Q, A is identically zero and A ± take values in the binary set 
{0, 1}. We refer to such special solutions as two-phase solutions. 

Assuming the smoothness of the interface separating the two stable phases, Evans and PortilheirollH 



took advantage of the entropy condition ( |10| ) to determine pointwise conditions to be satisfied on 
such transition interface by a piecewise smooth solution, in the same fashion as in the case of 
hyperbolic conservation laws. Precisely, given the smooth curve y := {(£(?), t)\t e [0, T]}, let 

/?* := {(x, f) ; te [0, T], ±(x - £(t)) > 0}, 

and assume u to be smooth outside y and such that, for any t e [0, T], 

u(R~) c (-00, b], u(R + ) c [a, 00), 
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and assume that there exist finite limits 

lim «(f(f)±i7,f)=: «(f (*)*,*), 



<9w 3w , 

lim — m^VJ) =■ -z-i&fj). 

r]->0 + OX OX 



Then condition ( |10| ) can be translated in a pointwise condition to be satisfied along the phase 
transition curve y. 

Proposition 2.1. [5] A piecewise smooth weak solution u with the above described structure sat- 
isfies the entropy condition ( |10[ ) if and only if: 

i. the function u is a classical solution of ([T]) outside the curve y; 

ii. along the curve y there hold 

d(p(u) 



(11) 



[<f>( U )]y = 0, [U]y + 



dx 



= 



(12) 



where [f] y := f(£(t) + , t) - f(£(t) , t) denotes the jump of f at y. 
iii. along the curve y the following entropy conditions are satisfied 

f(o>o if mm,t)=A, 
f(o<o if mm,t) = B, 

U'(0 = ifmm,t)e(A t B). 

Roughly speaking, a two-phase solution can face alternatively two different possible evolutions: 
the steady interface problem and the moving interface one. 

i. Steady interface. If (f>(u)(g(f), f) e (A,B), then the entropy condition ( fT2] l implies that the 
transition interface y is steady. Moreover, conditions ( fTT| ) become 

d(f>(u) 



dx 



= 0. 



In this case, the solution can be obtained by solving ([T]) at the left and at the right of y, with the 
condition of a C 1 -connection for cp(u) at y. 

ii. Moving interface. If cp(u)(^(t),t) e {A,B} then the entropy condition ( fT2| ) implies that the 
transition interface y is allowed to move, becoming a free interface of the problem, with direction 
of propagation depending on the value cf>(u)(^(t), t). To fix ideas, let <f>(u)(^(t), t) = A. The values 
of u at the left-hand and at the right-hand side of y are hence given by 

u{^{ty,t) = c, u(t(t) + ,t) = a. 



The second condition in ( |TT| ), to be read as a Rankine-Hugoniot condition, is needed to determine 
the speed of propagation of the free-interface. 

A result on local existence and uniqueness of such kind of solutions has been proved in Ref. |[T4ll . 
showing that the approach by Plotnikov is indeed sufficient to guarantee well-posedness when re- 
stricted to solutions with pure and stable phases. 
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Going back to the derivation of the complete model equation, it is natural to ask if the two-phase 
description for the limiting solutions to ([2]) as s — > + still holds when considering the limiting 
solutions to ([5]> as s, t — > + . We are not aware of any result in this direction and we regard the 
problem of determining an appropriate limiting dynamics for (|5]> as a very challenging direction 
to investigate. 



Riemann problems. Among other initial-value problems, the Riemann problem for ([T]), deter- 
mined by an initial data of the form 

if x < 0, 

w + x > 0, 



(13) K°(X) = 



with m* e R given, is particularly noteworthy because of its invariance with respect to the rescaling 
(x, t) \-> {Ax, A 2 t). As a consequence, for any € R, the solution of the Cauchy problem ([T])-([T3|) 
has the form 

(14) u(x,t) = f(Q where f := — 

where the function / is such that f(±oo) = u ± . Discontinuity curves in the (x, ?)-plane have the 
form x = | yft for some f € R. For £ # the function / is solution to the following (second order) 
ordinary differential equation 

(15) 0CO" + ^/'=O. 



At any jump point the following relations, obtained from conditions (|TT 
(16) [0(/)]f = O, + if [/]| = 0, 

have to hold. A detailed analysis of such self-similar solutions to ([T]) has been fulfilled in the case 
of general initial data (hence, also in the unstable region) ifTUll . Here, we focus on the case of a 
piecewise linear diffusion function <p 

(f>{u) - (f>~{u) :— m~u + q~ for u <b, 

(17) 

4>{u) - (p + {u) := m + u + q + for u> a, 
for some m ± > and q ± € R. In this case, it is possible to assemble an explicit formula for the 



solution satisfying the conditions described in Proposition 2. 1 Such solutions will be useful to 
test the algorithms proposed in the subsequent Sections. 

In order to give an explicit formula for the solution to the Riemann problem, let us introduce 
the functions, related to the fundamental solution of the linear diffusion equation in the half-space: 
given m > 0, we set 



1 r*^ 1 f 

E%£) := -= e-y 2 " m dy, E + J£) := — = 

y4nm J -co V4 n m J£ 
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The functions F* satisfy the properties 

F±(±oo) = 0, = 1 and £-+£+ = 1. 

Without loss of generality, we consider the case u~ e (-00, b] and w + e [a, +00). 

Proposition 2.2. . Let (p be given as in (T7l. 77ze Riemann problem for ([T]), determined by the 
initial datum (|13[), /jas a two-phase solution with a steady interface if and only if 



: (f>(u ) + 



(18) 



yrrr + ym ynr + ym 

In this case, the solution has the form (fT4]l where 



m ■= 



{<p-T\g(j;)) 



and 



(19) 



gig) ■-- 



<P(u-)E + (Z) + <P(u + )E-(0 - 



(P(u)E + m+ {Z) + <t>{u + )E- m+ {i;) 



Jm^ 



4>(u + ) € [A, B], 



Jnr + ym 



/m T - ym 



^>0. 



ynv + "vm 

Proof. Let us look for the self-similar solution / in the form 

u+C-E- m ^) f<l 
u + -C + E + m ^) f>l 
where C ± have to be determined. In terms of the function 

>(«-) + K~ E~{0 €<l 



g(g) := = 



<f>{u + ) - K + EU{) 



Z>1 



where k* := m* C*, the transmission conditions ( fT6| ) give a linear system for the unknowns 

'£;L(D* + + jr-(D*- = M 



-| 2 /4m + 



m + 



K + + 



y^l 



m 



Setting 



A(|) :=F+ (f) 



,-eiAm- 



lm 



+ 



m 



-f 2 /4m+ 



k = y^r^[u] 



nv 



E + m+ iO 



-|| 2 (l/m + +l/m-) 



(m + - m ) - 
F + (& + F_(0 + \ _ f F + (0 F_(0 , 



where F,-(£) := yTn' e^ 2/4m ' E 1 .(£), with / e {-, +}, the values of k* are given by 



k + = 



[0] + ^V(« + ) - ^(« + )) 



A(|) 
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e -f/4m + 



In the case E, = 0, there hold 

A(|) 



m + A(£) 



[0] - ^(<f> + (u~) - *"(«")) 



2 



so that the function g is given by (f!9|). Such a formula defines an entropy solution of the problem 



if and only if condition iii. in Proposition 2.1 is satisfied, that is, 



in 



1 _ 1 . 1 Vm + 

-4>(u ) + -<p(u + ) - -—= = 

ILL yj m + + yj m 



that is equivalent to ( 18 1. 



[<t>] e [A,B], 



If ( fT8| ) does not hold, the Riemann problem exhibits a moving interface. For the sake of sim- 
plicity, let us restrict the attention to the case m + = trT - m. Then the values of at* reduce to 



= [0] : 

and the function g can be rewritten as 



(20) 



gig) = \<Ku) + -^:(<f - q + )tet /4m E-(i)x [l+ JZ)\E+(f) 



where x, denotes the characteristic function of the set /. Condition ( [T8] ) corresponds to the request 

l -{<p{u) + cp{u + )) e [A,B]. 

If \(<fi(u~) + <f>(u + )) > B (or < A), the relation g(£) = B (= A, respectively) is an implicit relation 
for the value £ and the transition curve for the corresponding solution to the Riemann problem is 
given by x - £(t) = | y/i. 

3. Semi-discrete approximation of the pseudo-parabolic equation 

Next, we aim at designing appropriate numerical schemes for both equation ([TJ), considered in 
the setting described in Section [2j and the pseudo-parabolic equation ([2]). First of all, we concen- 
trate on the latter, considered in a finite interval / c R, with homogeneous Neumann boundary 
conditions. ifTTTl 

du 



(21) 



dx 



xedl 



= 0. 



Different types of approximations of pseudo-parabolic equations have been considered in the liter- 
ature (finite-element methods Reff. JTl [T3j [30j [8j [7]| , spectral methods Ref. [23], finite differences 
Ref. IfTTTl ...). Here, we consider a standard finite-difference space discretization for equation Q, 
that appears to be continuous at e = 0, so that the limit £ — > + can be analyzed by simply setting 
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e = directly in the algorithm. Formally, this amounts to an exchange of order of the two Vanish- 
ed 

ing limits h — > and s — > . Because of the presence of the ratio — , no uniformity of one limit 
with respect to the other parameter into play may hold, in general. Hence, there is no guarantee 
that such an approach leads to the correct limiting solution. 

The semi-discrete scheme. Let / := (0, 1). Given h :- 1/(7 - 1) > for some J > 5, let us set 

xj :- j 'h for j - I,,.., J. Let U — U(t) - (U\(t), . . . , U j(t)). Let us approximate the second order 
derivative with respect to x with 



d 2 u 
dx 1 



(Xj,t) 



u(xj + i,t) - 2u(xj, t) + u(xj-i,t) 
h 2 



+ 0(h 2 ). 



The boundary conditions ( |2Tj ) translate into the conditions 

ui(t) = u 2 (t), uAo = uj-m 

For the sake of simplicity, we denote abusively (<p(U\), . . . ,<p(U j)) T by (p(U). Equation ([2]) can 
thus be approximated by the ODE system 



dU 



(h 2 I + £A) — = ~Ad>(U), 
' at 



where I is the identity matrix and 



(22) 



' 1 


-1 


o 


. ... o 





' 


-1 


2 


-1 '•. •• 


. ... o 











-1 


2 "-. ' 


. ••• 














o ... • 


. '•. 2 


-1 











o 


• '•■ -1 


2 


-1 


I o 





o 


. ... o 


-1 


1 , 



The matrix A being an irreducible Hessenberg symmetric matrix, A is diagonalizable and all its 
eigenvalues are simple. By Gerschgorin-Hadamard's circle theorem, the spectrum <r(A) of the 
(symmetric) matrix A is contained in D(2, 2) := \X € C : \X — 2| < 2}. Let us denote the 
eigenvalues of A by 



(23) 



. . 2 ( U-l)n 

4 sin 

2/ 



je{l,...,J], 



as can be inferred from the computation of the eigenvectors. As a consequence, the matrix A is 
positive semi-definite and the matrix h 2 I + e A is invertible for any value of positive e and h. Note 
that A j converges to 4 as J — > oo, so that the upper bound given by Gerschgorin-Hadamard's circle 
theorem is optimal. 



14 



The semi-discrete scheme can be rewritten as 
(24) = - (h 2 I + eA)~A0(£/) =: -A h 2 E (f>(U). 

Remark 3.1. Since A h 2 s is a polynomial fraction of the matrix A, the spectral theorem implies at 
once that A and A h 2 e share the same spectral decomposition, that is their eigenvalues are simple 
and their eigenvectors are associated respectively with Aj and R(sAj/h 2 )/e, for any j e {1, . . . , /}, 
where R : x e R + i-> ;t/(l + jc). Since for any x > 0, < /?(X) < min(l,jc), the spectrum of 
A/,2 e is bounded from above by min(l/e,4//j 2 ). Thus the numerical effect of considering A^2 >e 
in place of h~ 2 A is that the large eigenvalues of A are filtered. Note also the difference becomes 
more apparent when dealing with data in the unstable region, cp' < 0, since for such values the 



semi-discrete scheme ( 24 1 exhibits exponentially increasing behavior at the linearized level. 



Following Remark 3.1 and the fact that A h 2 E has a regular limit as e — > + , we now focus on 
the limit of the semi-discretized system 

dU 

(25) — = -AftU), 

CIT 



where t = t/h . In the original variable t, system d25j) is stiff as h — > . Note that, since -A o <p 



is autonomous and globally Lipschitz continuous, the Cauchy problem for system ( |25) with an 
initial condition admits a unique C 1 solution defined globally. The study of the behavior of this 
semi-discrete solution requires the knowledge of both L°° and L 2 properties of System (25 1. Since 
(p is assumed to be piecewise linear, this system can be seen as 

(26) ^-(t) = -M(t)U(t) + W(t), 

CIT 

where M(t) is a J x / matrix and W(t) is a vector of size J, M and W being piecewise constant in 
time: the jumps occur when a point changes phases, that is when the interface moves. 

Remark 3.2. From now on, we will concentrate on a particular symmetric piecewise linear <p given 
in (|7]) such that ^(s) = 2s + 3,a = -b= 1, d = —c = 2 and B - -A - 1 (see Figure|2]). In this 
case, condition ( p~8] > becomes u~ + u + € [-1, 1]. The numerical analysis we will perform makes use 
of the remarkable spectral properties of the matrix A. If cp has different slopes m ± > on the left- 
and right- hand sides, then one can apply the following analysis keeping in mind that A should be 
replaced by AD, where D is a diagonal matrix of the form diag(m~, . . . , m~, m + , . . . , m + ), so that 
A and AD are similar matrices and share the same spectral properties. 



Riemann initial data with steady interface. To analyze the semi-discrete algorithm ( |24| ), we 
consider the Riemann initial data 

u~ x < 1/2, 

u + x> 1/2, 

where e <S ± are chosen such that u~ < b < a < u + satisfy ( fT8j ), so that the interface should not 
move and both sides should remain in their respective stable phases. 



u°(x) 
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First of all, let L e {1, . . . , J} be such that the initial datum U° = (U^, ...,U°j) for the system 



( |25j ) satisfies 

(27) U°<b<a<U° h Vj=l,...,L, h = L+\,...,J, 

The behavior of the solution £/ can be easily determined since </) is supposed to be piecewise linear 
and symmetric, as given in Q. Indeed, the function V = (p(U) = 2U + (3, ... ,3, -3, . . . , -3) T 
satisfies, at least locally in time, the system 

dV 

(28) — = -A V, 

as 

where s = It - 2t/h 2 , so that, as long as the components of the solution U are localized as in 



( |27j ), there holds 

V(s) = exp(-sA)V(0). 

Since exp(-sA) = lim (I + sA/n)~ n , ||(I + sA/n)' 1 ]^ = 1 and (I + sA/n)~ l e (R + ) 7x/ , we infer 



n— >+oo 



|| exp(-sA)||oo = 1 and exp(-sA) € (R + ) 7x/ . 

Consequently, if V y (0) € [-1,1] for any j, then Vj(s) € [-1, 1] for all j € {1, . . . , J} and for any 
s > 0. If Vj(0) <f. [-1, 1], that is u + > 2 or iC < -2, and u + + u~ e [0, 1] (resp. u + + iC e [-1,0]), 
one can get back to the previous case by dividing U° by u~ (resp. u + ). 

Hence, there is no change of phases and M(t) = -A for all t > 0. This behavior is of course 
compatible with condition ( fT8] >. 

Let us now study the asymptotic behavior of the solution U. The value belongs to o~(A) and 

Ar - where r\ = (1, . . . , l) r / V7. Therefore, the matrix Pi := r\r T x - 1/J, where 1 is the 

/ x / matrix composed of ones, is the eigenprojection corresponding to the eigenvalue 0. Hence, 

1 J 

the matrix A admits a spectral decomposition of the form A = — I- 2 At Pj, where A; are the 

J ;=2 

eigenvalues of A and Pj are the corresponding eigenprojections. Therefore 

1 J 

(29) e~ Al = - + Y i e~ XjS Pj 

and a solution to ( |2"8) ) with initial data V° converges exponentially fast to 

PiV° = (- T f. Vj)(l....,l) r . 

;=i 

This means that the solution U converges to a Riemann-shaped steady state as t — > +oo. The 
solution to system ( [28] ) with initial data V° is given by 

V(r) = V° = (- } j] Vj) (1, . . • , If + 2 e-^Pj V°. 

7=1 7=2 
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If U° is a Riemann-type initial data with values (u , u + ), that is 



t/J = u 



Vj=l,...,L, U)l = u + Vh = L+l,...,J, 

the jump being at the middle point of the grid, we have 

<p(u + ) + <t>(u~) 



(30) V(t/hx z ) = 

Hence, there holds for all j € {!,...,/} 



(31) 



r t//v/i,2n 0(w + ) + <^(" ) 
hm Vj(t/hr) = = v 

/->+co J 2 



that is the value towards which cp converges at the transition interface £. Thus the solution U 
converges to a two-phase state 



(32) 



(<p-)-\v°°) if j<L, 



(</> + r\v°°) 



if j>L+l. 



Thanks to the spectral properties ( |23] > and ( |30] >, we know that this convergence is exponential of 
rate 



2A 2 = — sin 2 I — 



In 1 . 



Behavior at a transition point. Now, let us focus on a model situation in which one point is on 
the verge of the unstable phase 14, that is there exists L e {2, 1} such that the initial datum 
is given by 

U° - {b,...,b,d,d + 6) T , 
where the definitions of b, d are described in Figure [2j the value d is assumed at the (L + l)-th 
position of the vector V° and 6 > is a positive perturbation. The borderline point V° is attracted 
to the stable phase S + because of the jump conditions, but it has to cross the unstable phase 14 
first. Let us determine the expression of M(t) at least in a neighborhood of r = 0. We already 
know that the solution U is of class C . At first, note that U j(t) is in the stable phase <S + for 
j > L + 1 and that 

(25 



dr 



(0) = 0(t/,_i(O)) - 20(t/ y (O)) + <p(U j+1 (0)) = 



-26 




if j = L+ 1, 
if j - L + 2, 
otherwise, 



so that 



-l+o(r) Hj<L, 

2 + 26t + o(t) if j = L+ I, 

2 + 5-25t + o(t) ifj = L + 2, 

2 + 5 + o(t) ifj>L + 3. 

Consequently, one has (dUi/dTXr) = 45t + o(r), and therefore Ul increases and changes phases 
S- 14. Meanwhile, for j < L - 1, (dUj/dT)(T) = o(r), so that Uj(t) = -1 + o(r 2 ) and 
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u + s, 



(dUL-i/dT)(T) = -2St 2 + o(t 2 ). By recursion, Uj stays locally in the stable phase S~, for j < L-l 
and we can write the ODE system p5j ) linearly, at least locally in time in a neighborhood of t = 0, 

as 

(34) — = 

CIT 

withS := (0,. ..,0, -3,0,3,0,. ..,0) r and 

( 1 -1 

-1 2 -1 '-. ... . 
-1 2 '■. '•. 



:= 2 



















o ^ 





-1 2 -1 '■■ 

-. -1 2 1/2 *•. 

'•• -1 -1 -1 '•■ 

'■■ 1/2 2 -1 '•■ 

'■■ -1 2 -1 







2-10 


















-1 





-1 



■1 1 



where the "defect" is located at the (L - l)-th, L-th and (L + l)-th lines. The above expression 
allows us to give a formula for U(t), for as long as Ul stays in the unstable phase K, that is 
re [0,f]: 



(35) 



U(t) = exp(-rB){7(0) + exp(-sB)ds S. 



f' 

Jo 



In order to fully exploit ( |35j ) and predict the time f at which Ul shall leave the unstable phase H, 
let us describe the spectrum of B. 

Proposition 3.3. The matrix B is diagonalizable in R with J distinct eigenvalues which can be 
ordered as -1 < //-i < /jq = < fi\ < . . . < /j.j-2 < 8- In particular, B has only one negative 
eigenvalue. 

Remark 3.4. The negative eigenvalue of B is the one that allows the borderline point to go through 
the unstable phase II. 

Proof. The matrix B r can be written as B T = DA where D = diag(2, . . . , 2, -1, 2, . . . , 2). Let us 
prove that B r is diagonalizable. As mentioned before, A is symmetric non-negative and can 
be diagonalized in an orthonormal basis (ro, . . . , o_i). Let A := diag(0, A\, . . . , A := 
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diag(l,/ii, ...,Aj-i), which is invertible, and W be the unitary matrix the columns of which are 
(ro, . . . , rj-i) in that order, so that 



I 1 T \ 

w 



W = 



V7 



w 



V7 

with w e R 7 " 1 and A = WAW T where we denote by M the (1, l)-submatrix of size (7- 1) x (7- 1) 
of a 7 x 7 matrix M. Let us note that 

(36) jU-ij-i +WW T = ww T + W T W = Idj-u-u 

so that, in particular, W is invertible and W~ l = - - - W r . The matrix C := W r B r W = 
W r DlVW r AW = W T DIVA is similar to B T 



Let us now multiply on the left- and right-hand sides of C by A 1 ^ 2 and A l/2 : 
A l/2 CA- l/2 = A l/2 W T DWA l/2 



c T 
0/_ u A 1 ^ 1 / 2 

where c € M. J ~ l and A := W /T DW+2wH' r which is symmetric. Let us show now that A is invertible 
and has exactly one negative and 7-2 positive eigenvalues, so that B r (and B) is similar to a 
diagonalizable matrix with one negative, one zero and 7-2 positive eigenvalues using Sylvester's 
law of inertia. 



Using relation ( |36| ), one can rewrite A as 

A = 2Id + W T {D - 2Id)W. 

The definition of D implies that 

D - 21 d = -3E LL 

where El-h-i '■= (^i=i-i y=z,-l)l<; ;</-i> <5 being the Kronecker symbol. The spectrum of A is thus 

~_i /-3 
2 with multiplicity 7-2 and associated eigenspace W (Span(e j) j+l-i) an d — associated 

with W~ l eL-i- This proves that B is diagonalizable with simple eigenvalues with the announced 
signs. 

Applying Gerschgorin-Hadamard's circle theorem to B r , we have 

Sp(B)cD(4,4)UD(-2,2), 

so that we get the upper bound of the spectrum quite easily. To prove the lower bound, we 
need to examine more closely the computation of the eigenvalues: assume (ju, V) is an eigen- 
value/eigenvector of B T . Solving B T V = fiV leads to solving two linear recurrences of order 
2 that have the same characteristic polynomial with compatibility conditions. The computation 
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gives the following result: if /i € [-4, 0) U (0, 8) is an eigenvalue of B T , then 

X (M) -=-2(p + 2) [v J + - 2 (p) + v J S 2 (ju)} + (4 + rf2 - fi)) [v^fci) + v ] J V)] 
+ 3^[v^-V) + v^-V)] 

-0 

where v ± (ji) = 1 - ^ ± ~ i)- Note that > for a11 fe 6 M ' 

^>+^TU)('-r( 



2 



1 - 



p=0 

is a //-polynomial of degree k, so that ^ is also a /i-polynomial, of degree J + 1 . The possible 
eigenvalue // — 8 is not valid in^- since it corresponds to a double root in v ± . However, ,^(8) = 0, 
but a further computation shows that 8 is not an eigenvalue of B: it implies that the characteristic 
polynomial of B T is^- := x/($ - //). We already know that there is only one negative eigenvalue 
of x- Proving that^(-l)£(-oo) > will show that this negative eigenvalue is necessarily larger 
than -1. Since v±(-l) - 2* 1 , we get^(-l) - -(2 W + 2 J ~ 2L ~ [ + 2 2L+1 " / )/3 < 0. Thanks to the 
expression of v+ _1 (/i) + vf _1 (/i), one sees that it is positive as yu tends to -oo so that^(-oo) < 0. 
This ends the proof of Proposition [331 □ 



Since B is diagonalizable with simple eigenvalues, one can compute explicitly the spectral 
projectors: let fj. be an eigenvalue of B and (v, w) € (R J \ {0}) 2 such that Bv = /uv, B T w = /uw. 
Since w T v + 0, one can assume w T v = 1, so that the spectral projector associated to ji reads 
Pu = vw T /(w T v). The direct consequence is the following expression for U(t), t e [0, f], from 



(351: 



(37) U(t) = P (U(0) + tS)+ }^ Pje- Tti U(0) + S\ 

^eSp(B)\{0) ^ ^ ' 

so that, identifying the asymptotically larger terms and noting that, since Po = J D~ l ee T /(J - 1), 
there holds PoS = 0, we infer 

(38) U(T) = e- T ^P u Ju(0)-—s) + P U(0)+ V -P^S + o{e^ T ). 

V ' (U eSp(B)\{0} ^ 

An approximation of the exit time f is thus given by 



(39) 



log 


i [Pu8](L)\ 
a [P f/(0)](D Z /i€ Sp(B)\{0} \ ) 


log( 









Note that formulas (|33[)-(|37[)-(|38|>-(|39]> hold for general piecewise-linear cf>, changing B to a similar 
matrix and modifying S. 
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In conclusion, within the time interval [0,f], there is only one component, namely Ul, that 
lies in the unstable phase H and, as soon as it has entered the stable phase S + , the matrix of 



the dynamical system is again as in ( [28] >. It means that the interface between S~ and S + has 
moved to the left, as predicted by ( fT8j ). One is then faced to two independent constant coefficient 
heat equations, that can be treated as in the steady configuration case. If, on the right-side, the 
perturbation 5 is too large, Ul+i will at some point be larger than B and the situation will revert 
to that with a point in the unstable phase H. The system will then move the interface to the right- 
hand side again by making this point go through the unstable phase 1i, until the border is reached 
or the discretization is such that (P(Ul+i) is less that B. 

In the general case of initial data of Riemann type, satisfying ( |27] >, and such that \u~ + u + \ > 1, 
the semi-discrete differential system will make Ul converge exponentially fast to the "closest" 
border predicted by the entropy conditions ( fT2| ): if u~ + u + < -1, then (£/z,(f), UL+i(t)) goes to 
(c, a) and otherwise (Ui(t), UL+\(t) goes to (b,d). The rate of convergence to this intermediate 
state is again 2A2 = 8 sm 2 (n/(2J)/h 2 ~ In 1 . Once the border is reached, we are faced with a 

./— >+oo 



situation where the dynamical system is the same as in p4| ). The analysis previously developed 
still holds: the borderline point (Ul if u~ +u + > 1 or Ul+i otherwise) is the only one going through 
the unstable phase 14. Once it has reached the opposite stable phase, the following point (Ul-i or 
Ul+i), that has been monotonically changing (decreasing or increasing) goes to the border, then 



enters tl and so on, so that the interface moves as predicted by Proposition 2.1 in the continuous 
setting. 

4. Numerical experiments in the limiting regime e = 

Now that we have shown that the semi-discrete scheme plays the role of simulating either the 
convergence to a Riemann-type steady state or the movement of the interface, we turn our attention 
to performing concrete numerical experiments, by adding an appropriate discretization for the time 
derivative. 

There are two ways of fully discretizing the problem, remembering that equation ([T]) is fully 
non-linear and <p is not differentiable, so that Newton's algorithm cannot be used in a fully implicit 
scheme: 

- either discretize explicitly in time and space, that is introducing a time-discretization f = nAt 
with step At and using an explicit Euler scheme 

(40) U n+l = U n -^A<p(U n ), 
where (U^j^i j}^; 

- or use the previous analysis and discretize ((26]) with an implicit scheme, globally in time if the 
initial data is of steady Riemann type or else by estimating the exit/entrance time 

(41) = U n - ^(M(t n+l ) U n+1 - W(t n+l )) . 
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Using ( |40| ) allows to obtain a numerical solution without having to consider exit/entrance times, 
but the CFL condition for stability is 4A?//z 2 < 1 and the order of approximation is 1 in time, 
whereas ( |4T] > avoids L 2 and L°° stability conditions and increase the order of the approximation 
but requires the knowledge of M(t). Here, we implement both approaches and compare the results. 

First of all, we analyze the case of initial data in stable phases, in order to validate the schemes 



by comparing to the exact solution for the Riemann problem stated in Proposition 2.2 with the 
left-hand state u~ belonging to S~ and the right-hand state u + to «S + , with S ± defined in ([6]). 
Specifically, we choose the case (u~,u + ) = (-1, 1), for which the predicted behavior is that of a 
steady interface and continuity and differentiability of cp(u) at x - 0.5. The convergence error to 



mplicll scheme - u al time fl.005 - 2" (-11) h. 1/255 




1.5 2.0 2.5 



(a) u as a function of x (b) (p(u) as a function of u 

Figure 3. Numerical solution at time T = 0.05 with Riemann initial data (-1, 1). 



the exact solution given by ([20} can be expressed as 
(42) 



E 2 := 



A h Y(U". - u{t'\xj)) 2 = 0(At a ) + 0(hP). 



In the case of standard diffusion equations with regular initial data, e.g. taking the exact solution 
at small t for (u~, u + ) in the same stable phase, under the required CFL condition At = 0(h 2 ), the 
order of the explicit scheme is a = 1, ft = 2. The implicit scheme requires no stability condition 
and is also of order a = \,f} - 2, the counterpart being that we should solve a system at each step. 
Here, since the initial data are step functions, one cannot hope for such good estimates. 

Numerically, we show here that the order of approximation is ft = 1 in h for both schemes and 
between 0.5 and 1 in At for the implicit scheme. The higher-order region in At could be explained 
by the fact that, h being fixed, as we refine At, the solution is smoothed out faster. In conclusion, 
since in the steady case the system to be inverted is linear, tridiagonal, the implicit scheme should 
be preferred (see Table 1). 
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implicit (var At) 


implicit (var h) 




Explicit (At - 


= h 2 /4) 
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J 


error 


cpu 
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error 


cpu 
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error 


cpu 
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0.0043659 


0.02 




6. 


0.0264534 


1.76 


6 


0.0373874 


0.02 


14 




0.0034607 


0.06 




7. 


0.0132975 


1.94 


7 


0.0188008 


0.13 


18 




0.0033399 


0.62 




8. 


0.0066676 


2.31 
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0.0094283 


0.16 


20 


9 


0.0033395 


3.05 


20 


9. 


0.0033395 


3.05 


9 


0.0047216 


0.80 


22 




0.0033394 


9.98 




10. 


0.0016730 


4.48 


10 


0.0023638 


3.87 












11. 


0.0008409 


7.71 


11 


0.0011850 


25.54 












12. 


0.0004287 


14.37 


12 


0.0005981 


193.12 












13. 


0.0002299 


33.90 


13 


0.0003099 


1971.36 


















14 


0.0001749 


17407.04 



Table 1 . Comparison of L 2 errors of <p(u) and execution times (cpu) for the implicit and 
Explicit schemes wit At = 2~ K , h = 2~ J . 



Next, we consider a Riemann-type initial data with values (u~, u + ) = (-2, 3). Since <f>(-2) = -1 
and 0(3) = 3, the function <p(u) is discontinuous at time t = with a jump located at the middle 
point x = 0.5. In short time, as expected, the function <p(u) is regularized and, at x = 0.5, takes 
a value in the interval [A, B] = [—1,1] depending on the choice of u~,u + (see Figure [4]). As a 
consequence of such a regularization of <p, the profile of u is forced to bend at the left-hand and 
right-hand side of the phase transition point. From now on, diffusion in the two stable phases will 
flatten the profile of u on the two phases and, at the same time, the interaction between the two 
stable phases will determine the position of the transition interface and of the corresponding value 
of <p(u) on the interface itself. 



Implicit scheme - u al time 0.005 -dt^E' (-11) ■ h = 1/2S5 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.0 0.9 1.0 -3 -2 -10 1 2 3 4 



(a) u as a function of x (b) (p(u) as a function of u 

FIGURE 4. Short time evolution of the Riemann data (-2, 3). The two graphs show the profile of 
u and <f>(u) at the final time t = 5 x 10~ 3 . 

For large time, the solution u converges to a Riemann-shaped steady state and the profile of (p(u) 
tends to a constant exponentially fast, as shown in Figure [5] As the continuous theory predicted 
(p~8]), since cp{iC) + cp{u + ) - 2 e [2A, 2B], the interface did not move. The same asymptotic 
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behavior is shown by different choices of initial values (u , u + ). As long as we were dealing with 
steady cases, we were able to use the implicit scheme without further study. 




(a) uasa function of x (b) error as a function of x 

FIGURE 5. Large time evolution of the Riemann data (-2, 3). Letting the time go as far as t = 0.5 
ensures an error of order 1(T 4 according to {33}. 

A higher level in the jump between the initial states u~ and u + leads to a movement of the 
transition interface separating the stable phases, as shown in Figure [6} corresponding to the initial 
data (-2, 4). We used an explicit scheme to treat this case since a point might have to go through 
the unstable phase. 

The analysis of the movement of the transition interface £ starting at x = 0.5 at initial time t = 
is very delicate. First of all, we have to define the position of the interface. Here, we consider 
the interface £ in the semi-discretized scheme to be determined by finding the maximal (resp. 
minimal) value of j such that Uj < — 1 (resp. Uj > 1). Here, the position £ = £(t) of the interface 
is given by £(t) = k{t) ■ h, k(t) being the index localizing the left-hand side of the interface. The 
explicit scheme allows us to get an approximation of the movement of the interface (see Figure[7]), 
from which one recognizes the parabolic shape predicted by ( fT4"l ). However, the use of this scheme 
yields oscillations. 

The space derivative 4>{u) x of the function cf>(u) varies very rapidly close to the interface <f, 
consistently with the presence of a movement of the interface. Nevertheless, the algorithm finds 
as value 0(m)|^ an uncorrected level (in the limit e — > + , left-moving interface are allowed if 
and only if <p(u) = B = 1). In fact, in the continuous limiting model the unstable phase (a, b) is 
prohibited and the solution u jumps instantaneously from one phase to the other. On the contrary, 
in the semi-discretized algorithm p5| ), the transition from one phase to the other is accomplished 
by smoothly passing from values under the level a = -1 to values above the level 6=1, thus 
passing all way through the spinodal region {a, b) (see Figure [6jc)-(d)). Correspondingly, the 
function <p decreases its value from B to A and then increases again to B. Such a transition is faster 
as h — > + , but it is always (incorrectly) present. 



(a) u as a function of x 



(b) <f>(u) as a function of u 




(c) <f>(u) as a function of x (d) zoom of around the interface 

FIGURE 6. Small time evolution of the Riemann data (-2, 4). On the right-hand side, the graph in 
the phase-space (u, <f>(u)) with corresponding values of the solution, where a point in the spinodal 
region is clearly seen. 




(a) Movement of the interface (b) Values of —\d x (<p(u)y\l\u\ 

Figure 7. Riemann initial data uT = -2, u + = 4 with h = 0.01. 
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5. Two-phase scheme 

In this Section, we consider the problem of determining a numerical scheme for the limiting 
equation ([T} together with the entropy conditions described in Proposition 2.1 The algorithm is 



based on the idea that initial data taking values only in the stable regions give rise to solutions with 
values in the stable regions matched in order to satisfy the appropriate transmission conditions. In 
the interior of any region, where the solution is in one of the stable phases, a diffusion equation 
has to be solved; and the procedure has to be completed with an appropriate description of the 
transition interfaces behavior, based on the entropy conditions. 

At each step, additionally to the values U" e R J giving the approximate value of the solution, 
we also track the non-integer position (" of the transition interface determining the passage from 
one phase to the other. The index of such position is denoted by j". 

Description of the algorithm. Let 0* : R — > R be two strictly increasing extensions of restric- 
tions fe s and (f>\ respectively (where, for the sake of simplicity, <p is piecewise linear and 
symmetric). For j £ {2, 1}, let us introduce the projection matrix IL; 



);•_! 

l 2 
7 _ 7 -_i 



where Oy and lj denotes the null and the identity j x j matrices, respectively. Finally, given A < B, 
let us introduce the truncation function T defined by 

(43) T(U) := max{A,min{£/,fi}). 

The algorithm reads as follows. Let U € R J and j, e {1, . . . , J} be the values for the approximate 
solution and approximate interface location. Let C := C(U,j*) be the transition value suggested 
by ( pi) , i.e. set 

C = C(UJ,) := l - + <P + {Uj, + 2>) ■ 

If C £ [A, B], that is the transition value is not in the hysteresis interval, we use the truncation 



function ( |43| > to transform it in an admissible value. Specifically, if C > B, then T(C) = B; if 
C < A, then T(C) = A. Correspondingly, the values of U at j t and j* + 1 are mapped in the values 
((f>-y l (T(C)) and (4> + y\T(Q), respectively (i.e. b and d if C > B, and c and a if C < A). If 
C e [A, B], the transition value is in the hysteresis loop and the new values for U at j\ and /» + 1 
have to be determined, imposing the transmission conditions 

[000] = [<p{u) x ] = 0, 

that is, in discretized form, 

<KUj.) = tfUj. + i), <p{Uj,) - 4>{Uj,- { ) = cf>(U jt+2 ) - 
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These conditions give 

(44) tfUj.) = <P(Uj, +1 ) = X - (tfUj.-i) + <KUj. + 2)) = C(U,j.). 

Hence, let us set 



¥(U, j.) := (0, . . . , , (c/>-)-\T(C)), (cf> + T l (T(C)), 0, . . . , 0) 



(45) 



j-th element 

Then, setting U n :=¥{U n , ft) + (I - IF*) U'\ the scheme can be written as 

(n n u n+l =F(u n ,ft), 



(46) 



(l-U fl )U n+l = (I-U fl )\U n 



At 



A(f>{U n ) 



where A is the matrix introduced p2| ). Note that the scheme ( |46| ) coincides with ( |40| , e = 0, apart 
from the definition at the transition points j* and j* + I. 

To define the transition interface sequence, given U e R J and j* € {1, . . . , /}, let us define the 
values for the (discrete) time-derivative of the interface £ as follows 



m = 



i 



(rr 1 (T(Q)-(rr 1 (nQ) 



<KUj.+2)-T(C) T(C)-<KUj.-ir 
h h , 



that is 



^(C) = 



2(T(C) - C) 



((P)-\T(C))-(cf>-)-\T(C)))h 

If C e [A, B], then T(Q = C, and f = 0; if C > B, then T(C) = B and £ < 0; finally, if C < A, 
then T(C) = A and £ > 0. Thus, the entropy jump conditions are satisfied. The transition interface 
algorithm is defined by 



(47) 



and 



•n+l 
7* 



+ 1. 



Next, let us fix our attention on the scheme ( |46| ) in the steady interface case (i.e. C(U n , ft) e [A, B] 
for any n). In this situation, it is possible to consider a corresponding semi-discrete scheme. Setting 



V = <p(U), relations ( |44| ) becomes 

(48) Vj,=Vj, +1 = ^(Vj.. 1 + Vj. + 2), 

so that the original system for the variable V = (Vf, . . . , Vj) can be reduced to a system for 

V :- (V\, . . . , Vj m -i, Vj,+z, . . . , V J ). The system to be satisfied by V, taking into account ( [48] ), is 
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where s = 2t = 2t/h 2 and A is the symmetric (J - 2) x (J - 2) matrix defined by 



(50) 



1 

-1 




-1 





-1 

2 





















2 -1 

-1 3/2 -1/2 '-. 

-. -1/2 3/2 -1 

-1 2 







-1 


















-1 





2 

-1 



-1 
1 



(here the lines (. . . , 0, -1, 3/2, -1/2, 0, . . . ) and (. . . , 0, -1/2, 3/2, -1, 0, . . . ) correspond to the in- 
dices j t - 1 and j t + 2). The matrix A defined in ( |50"| has spectral properties similar to the 
ones of the matrix A. Specifically, A is semi-definite positive, e <r(A) and An = where 
h = (1,. • • , l) r / V/-2. The matrix P Y := h ® h = 1/(7 - 2), where 1 is the (J - 2) x (7 - 2) 
matrix composed by all ones, is the eigenprojection corresponding to the eigenvalue 0. As a con- 
sequence, solution to ( |49l ) with initial data V° converges exponentially fast to 

J-2 



that is the solution U converges to a Riemann-shaped steady state as t — > +oo as described in ( 32 1 



Numerical experiments for the two-phase scheme. Next, we perform numerical experiments 



for the scheme ( |46| ) with the same Riemann-type initial data considered for the scheme ( |40| ) specif- 
ically with (u~, u + ) = (-2, 3.5) and (u~, u + ) = (-2,4). Results are shown in Figure[8] 

In term of the solution profiles u (or <p(u)), the results obtained with the present approach and 
the one considered in the previous Section are consistent with each other. On the contrary, in term 
of position of the moving interface, the behavior is very different, being much more regular and 
non oscillatory in the case of the algorithm ( |46l > (compare Figures |7Jb) and [9]). 

Indeed, the two approaches differ strongly in the determination of the interface position. When 
dealing with the algorithm ( |40] >, the position of the transition interface is determined by the last 
index j such that u; < b. As a consequence, the movement of the interface is always concentrated 
at specific iterations and have an amount of multiples of the space-discretization size h. The 
algorithm proposed for the limit model ([T| is based on the assumption that both the solution u and 
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0.5 0.6 



0.0 0.1 0.2 0.3 0.4 0.5 



0.7 0.B 0.3 



(a) u = -2, u + = 3.5 



(b) u = -2, u + = 4 



FIGURE 8. Evolution of the Riemann data. The two graphs show the profile of u at time t = 0.2 
and h = 0.01. The choice of different symbols helps keeping track of the initial position of the 
interface. 




(a) Movement of the interface 



(b) Values of -[dAd>(u))]/[u] 



Figure 9. Riemann initial data u = -2, u + = 4 with h = 0.01. 



the interface position £ are unknowns of the problem. Formula ( |47j ) permits to define a smoother 
location of the transition interface since its position is determined by calculating the values at 



each step. The discrepancy between the two approaches is shown in Figure 10 for the Riemann 
problem with initial data (-2,4). The two-phase scheme and the explicit scheme are both of order 



1 in h (see Fig. Ill, but the error for the two-phase scheme (compared to the exact solution 



computed in ( [20] ), the constant ^ being computed with Newton's algorithm) is less (see Figure 
TO] ) and the interface £ is better approximated and smoother (see Figure fT). However, due to the 
intermediate computations, the numerical cost is also higher than the one of the explicit scheme 
(see Table 2). 



(a) Solution of explicit and two-phase schemes 



(b) Zoom 



Figure 10. Comparison of the numerical solutions of the explicit and two-phase 
schemes with the exact solution for (zT , m + ) = (-2,4) at t — 0.005 for h — 1 /255 



L2 error ol phi(u| in log scale asafjrction of di for dU h " 2/4 at T= 0.005 L2 error of phi(u) in log scale as a function of dx for dt- h " 2/4 at T- 0.005 




5 6 7 S 9 10 11 10 5 6 7 S 9 10 11 

-it>ga(*o -iog2{d3t) 



(a) Two-phase scheme: error as a function of the mesh (b) Explicit scheme: error as a function of the mesh 
size h h 



Figure 1 1. Convergence error for <p{u) in L 2 norm 



Two-phase scheme 


Explicit 


J 


error 


cpu 


J 


error 


cpu 


6 




0.028548 


0.23 


6 


0.0600245 


0.13 


7 




0.0133988 


0.34 


7 


0.0300282 


0.26 


8 




0.0065732 


1.16 


8 


0.0150125 


1.05 


9 




0.0032091 


7.18 


9 


0.0075057 


6.02 


10 




0.0015852 


64.73 


10 


0.0037539 


52.31 


11 




0.0007498 


919.69 


11 


0.0016060 


723.27 



Table 2. Comparison of L 2 errors of <p(u) and execution times (cpu) for the two-phase 
and explicit schemes for (u~ , u + ) = (-2, 4). 
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